Bethe Ansatz approach to the pairing fluctuations in the mesoscopic regime 
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We review the exact treatment of the pairing correlation functions in the canonical ensemble. The 
key for the calculations has been provided by relating the discrete BCS model to known integrable 
theories corresponding to the so called Gaudin magnets with suitable boundary terms. In the present 
case the correlation functions can be accessed beyond the formal level, allowing the description of 
the cross-over from few electrons to the thermodynamic limit. In particular, we summarize the 
results on the finite size scaling behavior of the canonical pairing clarifying some puzzles emerged 
in the past. Some recent developments and applications are outlined. 
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I. INTRODUCTION 

When a small attractive interaction is switched 
on in a Fermi gas, bound states are formed non- 
perturbatively fil . This is the essence of the BCS pairing 
phenomenon [2j-l5| , leading ultimately to a phase transi- 
tion as a result of a competition between the kinetic en- 
ergy and the tendency of Cooper pairs to condense 0, [|J . 
The transition is usually described as a gauge symme- 
try breaking where the phase of the order parameter ac- 
quires a fixed value. Nevertheless it should be observed 
that because phase and number are conjugate variables, 
a definite value of the phase can be consistently reached 
only in open systems, where the number of Cooper pairs 
can fluctuate. The thermodynamics of macroscopic sys- 
tems, however, is not affected by equilibrium fluctuations 
(grand-canonical and canonical statistical mechanics are 
equivalent in the thermodynamic limit). Accordingly, for 
macroscopic systems the BCS condensate is a quantum 
coherent state of Cooper pairs characterized by a well de- 
fined order parameter, which is the Cooper-pair binding 
energy A p], Q. Away from the thermodynamic limit, 
it becomes delicate when speaking of the existence of 
the BCS state if the BCS pairing energy overlaps only a 
few energy electronic levels (the electrons level spacing 
5e is inversely proportional to the volume of the system) . 
The point was famously remarked by Anderson with the 
question: "What is the size limit for a metallic particle to 
have superconducting properties?" 0- This conceptual 
challenge has been revived significantly by experiments 
on isolated metallic grains of nanoscopic size 0. The 
crucial aspect in the experiments is that the number of 
electrons inside the grain is fixed due to the typically very 
low capacitance of the sample. Therefore the standard 
superconducting order parameter (ejej) exactly vanishes 
(the mean field approximation in the grand canonical en- 
semble is inappropriate). For BCS theory, the quantity 
playing the role of the "order parameter" is the pairing 



correlation function u 



i j 



where i,j are 



quantum numbers labeling electronic energy levels. 

Canonical pairing fluctuations were first analyzed nu- 
merically [Sf in order to study the physics of metallic 
nano-grains[l(|. In such mesoscopic regimes the pair- 
ing phenomenon appears as a cross-over region domi- 
nated by superconducting fluctuations sized by the ratio 
Because the gauge symmetry cannot be bro- 
ken at finite sizes, the system is characterized by these 
correlations (instead of a local observable). Interestingly, 
such a physical regime is shared with other important 
physical situations, notably in nuclei (Tlj . and more re- 
cently in systems of confined degenerate alkali Fermi 
gases [TH, El- Here, we comment that for 6 Li at quan- 
tum degeneracy, the mesoscopic regime can be achieved 
with 5e ~ A ~ 10" 12 cV (corres pon ding to a confine- 
ment frequency of u ~ 1 kHz) fl2l.ll4j . In turn, because 
the BCS order parameter is coherent on the length scale 
£ ~ fci?/A, mesoscopic fluctuations in the atomic gas are 
important for small enough cloud size R ~ £ (l5j. 

Since the system in mesoscopic cross-over regimes is 
characterized by strong quantum fluctuations, its physi- 
cal behavior is very sensitive to the approximations em- 
ployed and therefore exact results play an important role. 
The BCS model was solved exactly in 1964 with the sem- 
inal contributions by Richardson and Sherman [l6j . The 
strategy they adopted is in the spirit very close to the co- 
ordinate Bethe ansatz: first they considered the Cooper 
pairs as effective bosonic particles; then they were able 
to incorporate the constraint coming from the actual 
fcrmionic statistics into the many-body wave function. 
In 1967 Gaudin [ItJ realized that the Richardson solu- 
tion can be obtained with a variation of an approach he 
had pursued to solve the so-called Gaudin magnet fl8j . 
The Richardson solution remained unnoticed by the con- 
densed matter community until the late ninety's, when it 
was re-discovered rendering an understanding of the low 
temperature physics of metallic grains [l(| and later on 
for various applications in nuclear physics [l9| and cold 
atoms 

Even with the knowledge of the exact solution for the 
spectrum of the system, the computation of the corre- 
lation function is a highly non-trivial task. In essence, 
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the complications arise because in the correlation func- 
tions (ip\Oioc\'4>) the eigenstates of the Hamiltonian are 
not easily expressed in terms of the 'natural' states the 
operator Oi oc acts on. Therefore an exceedingly compli- 
cated combinatorial problem arises involvin g (s ums of) 
scalar products between Bethe eingenstates [21(. Using 
the exact eigenstates of the Hamiltonian, the diagonal 
pairing correlation function uu was obtained by Richard- 
son although it was not evaluated explicitly [22|]. A key 
progress for accessing the correlation functions exactly 
came from the observation that the BCS model belongs 
to the class of models that can be studied by the pow- 
erful techniques developed within the Quantum Inverse 
Scattering Method (QISM) [HJ. SpecificahYthc BCS 
model is indeed a twisted Gaudin magnet [23l - l26l ] related 
to disordered six vertex models [24|, H3 • The first assault 
on the problem with QISM protocols was done in Ref. [28| 
where various correlation functions have been computed. 
A major sim plif ication of the formulas was achieved by 
Zhou et al. [26| by application of the so called Slavnov 
formula for the calculation of the scalar products between 
Bethe states [29| together with the 'solution of the in- 
verse problem' [3 (1 |3l| . Correlation functions have been 
expressed as the sum of certain determinants. The ulti- 
mate progress has been achieved by Faribault, Calabrese 
and Caux who further reduced the complexity involved in 
the calculations, by applying certain reduction formulas 
for the determinants [32j. 

The aim of this article is to review the path we sum- 
marized above. In Section II we intoduce the BCS model 
and highlight its integrability. In Section III, we com- 
ment on the derivation of correlation functions by means 
of generating functions. Section IV is devoted to determi- 
nant representations of correlation functions making use 
of the solution to the inverse problem together with the 
reduction formulas worked out in [32j . Section V is fo- 
cused on the pairing amplitute in the canonical ensemble. 
Section VI has a short view on the thermodynamic limit 
of the model. Section VII presents further ramifications. 



II. THE BCS MODEL 

The BCS Hamiltonian is 



i,j=l 



c ji c jf M 



3=1 



g is the pairing coupling constant; the quantum numbers 
j € {1...0} label the single particle energy levels Cj 
which arc doubly degenerate since a € {t, 1} labels elec- 
tron spin states; Cj_ a and n^ a := Cj a Cj a are annihilation 
and number operators, respectively. 

In the following we explain the connection of the BCS 
model with the sZ(2)-Gaudin model. For this goal we in- 
troduce the fundamental realization of su(2) ~ sl(2) in 



4,t c L'. s ! ■■= ( c W+4,; c ^ - !)/ 2 - Thc s/ ( 2 ) " low - 

est" weight module is generated by the vacuum vector 



|0)„^|0) J= 
est" weight (sj 
of interest here [35| ) . 

S j := ( S j) 2 + 1 ( S ? S i 

The bilinear combinations Sj~S~ and S~ can be ex 
pressed in terms of Casimir and Cartan operators 



Sj\0}j = Sj\0)j, where Sj is the "low- 
-1/2 for spin 1/2, which is the case 
). The quadratic Casimir operator is 
SrS+),S^\0) j = s j (8 J -l)\0) j . 



S 3 s 7 



(2) 



Thc key observation is that the constants of the motion 
of the model Hamiltonian (JXJ) can be obtained by the 
QISM. The method follows an 'inverse' procedure to ob- 
tain a Hamiltonian that is intcgrable by construction. 
The starting point is to find a couple of matrices, L and 
R satisfying the Yang-Baxter equation 

R(u - v)L{u) ® L(v) = L(v) <g) L(u)R(u - v) , (3) 

where u € C is the spectral parameter. For the present 
case, the relevant matrices can be obtained through 

R x {u-z;?]) = + f(u - z,rj)(T® X , (4) 

where <x is the vector of Pauli matrices, and f{x,rj) := 
2rj/(jj — 2x) depending on the arbitrary parameter rj € R. 
The R matrix corresponds to X = er and z = in ([3)1 
while the Lax matrix L is obtained as Rg 

Rx defines the column-to-column and row-to-row scat- 
tering matrices, respectively, of the two dimensional six 
vertex model with inhomogeneities e = {ei, . . . , eo} [2l| . 
The monodromy matrix 

T(u\e) := Lq,{u - en) ■ ■ ■ Lx(u - ei) (6) 

satisfies the Yang-Baxter equation 

R{u-v)T(u\e)®T{v\e) = T{v\e)®T(u\€)R(u-v) . (7) 

Thc twisted monodromy matrix 



is then 



T(u\e) 



f(u\e) 



expa(77)5>f T(u\e) 



A{u\e) B(u\e) 
C(u\e) D(u\e) 



(8) 



(9) 



It satisfies the Yang-Baxter equation as well due to 



R 



0. The transfer matrix is the 

trace (over the auxiliary 2x2 space) of the monodromy 
matrix 



terms of electron pairs S i 



si 



(S- 



t{u\e) := tr {0) f(u\e) = A{u\e) + D(u\e) 



(10) 
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The latter is a generating function of integrals of the 
motion of the theory because they commute at different 
values of spectral parameters: [t(u\e), t(v\e)] — 0. The 
expansion t(u\e) = Y^Lo ^"^(^k) generates a hierarchy 
of integrable systems since 



J2 [i b (u\e),tMe)}=0 



(11) 



a=b+c=0 



The sum is on ordered partitions of a including b V c = 0. 

The first non trivial terms of the transfer matrix t(u\e) 
are 



1^ = 21+2^ — 
^ — ' a — i 



j'=i 



where 



¥3 



(12) 



(13) 



and 5j := (Sf , S|, 5?) are spin vectors; Sf = 1/V2(5J± 
iSj)). The operators (fTB")) define the twisted Gaudin 
magnet, where [tj,T(] = holds. By these integrals of 
the motion, the BCS model becomes connected with the 
Gaudin Hamiltonians Ej as the Hamiltonian ([lj can be 
expressed in terms of Tj as 

q n 
H = 2e J T J + 9 3 E T J Tl + const - (I 4 ) 

3=1 3,1=1 



which is manifestly integrable |23j ] : 

[ir,^-] = o,j = i,...,n 



(15) 



The exact eigenstates of the BCS model [Tg, [l8|] are 
obtained first by diagonalizing t(u) = A(u\e) + D(u\e): 

i(u\e)M{ui}\e)) = A(«, K)k) (16) 
where the Bcthe vectors are 

AT 



W{«i})> = l[B(u a )\0) 



(17) 



Then the eigenstates of Tj [23] are obtained by the quasi- 
classical expansion of the Eq. (fT6|) with \ it>) = \ipo)+r]\tl>i). 
results to be 



N 



\^) = l[S + (e a )\0) 



(18) 



where 



S 



S: 



and the rapidities u a 
Richardson equations 



uq + 7/e Q with e Q given by the 



s(e a ) 



1 

2^ 



N 

E 



a = 1, 



,N . 



e/3 



(21) 



Finally the eigenvalues of H are obtained via f| 14|) : 

#|£)at = S|£)at with £ = Ea=r e "- Throughout the 
article we will consider the half filling case TV = $1/2, 
unless it is stated differently. 

We observe that the operators (fT9j) span the infi- 
nite dimensional Gaudin algebra £7[sZ(2)]. The lowest 
weight module of C?[sZ(2)] is generated by the vacuum 
10} := ®SLi|0> i: S-(u)\0) = , S z (u)\0) = S (u)\0) , 
wheren s(u) is the lowest weight of Q[sl(2)]. We ob- 
serve that the integrability of the BCS model can be 
obtained as an algebraic property of (?[sZ(2)]. In fact 
the mutual commutativity of Tj descends from the rela- 
tion between t(u) := Y^=i T j /( u ~~ ^ £ j) ( T j are residues 
of t(u) in u = 2ei) and invariants (trace and quantum 
determinant H3) of Q[sl{2)}: 



t(u) = c(u) + s [2] (u) 
where c(u) is a twisted Casimir operator 



c(u) :=—S z (u) + 
9 



(22) 



(23) 



S z (u)S z (u) + - (S+(u)S-(u) + S-(u)S+(u)) 



and 



J2] 



(u) :=Y jS j/{u-2sjf 



(24) 



The property [c(u),c(v)] = is the origin of the integra- 
bility of the BCS model. Therefore finding the spectrum 
of the BCS model means finding the representations of a 
twisted Gaudin algebra (labeled by its Casimir operator) . 

We mention that the Richardson equations ([21]) arc 
intimately related to the algebraic structure of £[sZ(2)] in 
that they act as constraints on the lowest weight s(e a ). 
Thus, the difference between the BCS and Gaudin model 
amounts to a different constraint imposed on the lowest 
weight vector of Q[sl(2)] which leads to different sets £, 
£' (£' is spanned by the solutions of (f2Tj) when g — > oo) . 
We will use this fact to extend the Sklyanin theorem [34[ 
to the BCS model. 

In the next sections we will be focusing on the following 
M-point charge and pairing correlation functions (CF) 

M 



S± ( U ) : =E^T^T< 5Z («) : =E^^-' ( 19 ) (£\St...Si I \T) = (£\ II(n ife , t + n^-l)/2|.F) (25) 



3=1 

n 

s{u) :=^2 / Sj/{u-2ej) 

3=1 



3=1 



(20) 



fe=l 



{£\SrSf\T) = u i j{£,T) 



(26) 
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The vectors 



generated by 



TV 



(£\ := (0\]JS-(e a ), 

a = l 

N 



are exact iV-pair eigenstates of (Q]) (see Eqs. (fT8f . ([191) ). 
Here, we observe that the evaluation of the CF's (like 
([2"5f . ([25])) proceeds along the action of local operators, 
say Oioc, onto the Bcthe states, the latter involving a col- 
lective reorganization of the vectors of the local Hilbcrt 
space Hi oc . Therefore, for any fixed Bethe root u, Oi oc 
has to be commuted with B(u) to finally act on the vac- 
uum |0). This gives rise to a problem of combinatorial 
nature, whose solution is a non trivial task. In the next 
section we will look explicitely tame the combinatorics 
for the special case of ([25]) . ([26]). 



III. GENERATING FUNCTION FOR CF 

In [HI Sklyanin suggested how the combinatoric com- 
plications involved in the calculation of the correla- 
tion functions can be overcome resorting the Generating 
Function (GF) technique. He applied it [27j to the sl(2) 
Gaudin model [l8j . The key role in his approach is played 
by a reordering making use of the Baker-Campbcll- 
Hausdorff (BCH) formula for elements of the SL{2) loop 
group associated with the Gaudin algebra £/[sZ(2)]. In 
this section we exploit the common algebraic root of the 
Gaudin and BCS models to extend the Sklyanin theorem 
to the BCS model. The GF we will be looking at, is 



C{£,U,T) := (£\ J] S z (h)\T) 



(27) 



hen 



where the sets £, J- C C \ £o are (in general distinct) 
sets of solutions of the Richardson equations ([2"T[) ; £q := 
{2ej,j = 1...0}; ft C C\(£UJU£ ). The order of the 
correlation is the cardinality of ft: |ft|; \£\ and are 
fixed by the number of pairs N . For instance, the one 
and two point CFs correspond to |ft| = 1 and |ft| = 2 
respectively. 

Now we present the Sklyanin theorem for the GF of 
the si (2) Gaudin model and apply it to the BCS model. 
Therefore we need the notation of the set of coordinated 
partitions V = {Pi : I € l---\V\} of the sets £,J r ,'H 
(see Ref. [HI): the partition P G V is a set of triplets 
{Ti...T| P |}; the triplet T G P is T = {£ T , T T ,ftr) , 
where ^ T C £ J ^ Jt C J and ft T C ft such that 
\£t\ = \Ft\ > 0, \H T \ > 0. 

The GF has been evaluated for the sl(2) Gaudin model 
exploiting the BCH formula for the SL(2) loop group 



hen 

eeS 

where {S z (u), S ± (u)} € G[sl(2)} and <j>(x), T)(x), ip(x) are 
meromorphic functions for x £ C with residues <f>f , T)h , ip e 
respectively (U- This formula allows to rearrange the 
products between loop group elements in (|27p 



( ex P S 4>(x) ex P S kx) ex P 5 ^(x)> 

- < cx p s t( X) ex P s m ex P s t( X )> 



Sklyanin proved the following theorem [35 



Theorem. C{£ J 7 ) is given by the formula 

C{£,U,F) = {-l) N (28) 

E ( IT n T (\£ T \y nTl S(WTUHT)) J] < h ) 
■P \TeP ) he u P 



where 



S{C) = l/2m f s(z) Y[(z- y)~ l dz 
J r yec 



n T := -2\£ T \l(\£ T \ - 1)! , W T := &r U T T , 

andUp := 'H\U T( -p'Ht. C[£^ r hi^J : ^) is a polynomial in 
S with integer coefficients. 

Expression ([28]) depends only on the sets W := £UF and 
T-L [13, H3] ; for the Gaudin model W is a set of solutions 
of <[2T|) for g -> oo; for the BCS model W is a set of 
solutions of the Richardson equations pip for generic g. 
The scalar products of Bethe states (and their norms) 
are a corollary of the Sklyanin theorem (|28[) for Tl = 0: 
(^IJ 7 ) = C(£,0, J"). Its consent with the determinant 
formulas [3 Hi has been elucidated in Refs.jH H^- 
We point out that the GF ([2"T| has simple poles in the 
set £q. This will play a key role in the following. 

Correlation functions. The charge and pairing CFs 
are matrix elements of the su(2) Lie algebra (instead of 
elements of Q[sl(2)]) using vector states of Q[sl(2)]. The 
projection from the sl(2) loop algebra on its Lie algebra 
is performed by taking the residue of C(£,T-l, J 7 ) in the 
poles hi 
CFs d2S 



2e Jl for hi G "H, / £ {1...M}. The charge 



arc 



(£|^ . . . S* M \-T) = Jim (ft - f )C(f , ft, 7") (29) 
where H — > £q involve hi — > 2ej V/ and H — £q means 
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Ui(hi - 2e jt ). Using (28]) yields 



M 



{£\S{...S Z M \T) =(-l)^n s i 



(30) 



( II ri T S(W T ) J 



n 



n T |£ 2 



II (kr-lO 



where P t := {PeP: max|# T | = fc}; T fc := {T e P 
\H T \ = k}. The quantity S{W T ) is 
s(e) 



S(W T )= 



e£W T 



II (e-z) 



E 



*(d) 



y — 



n (<*-*) B ^ T (d-s/) n (d-x) 



(31) 
\ 



where e and d are elements appearing singly and doubly 
in Wt respectively. The pairing CF (|26|) can be extracted 
from C(£,% : F) where the vectors in (|27[) are (£| := 
(£\S-( Zl ) and |J) := S+Oa)!-? 7 )- Then m m (£, F) is 

u lm {£,F) = lim (zx - 2e ; )(z 2 - 2e ro )<7(£, 0,^(32) 

C(£,®,jF) is then calculated using the Sklyanin theorem. 
For I m formula (|32[) gives 



Ulm{£, J 7 ) = (-1) 



iV+1 



(33) 



where Z := {z\, Z2} and Zt is one of Z and m; 

"Pfc :={PeP: max|Z T | = k} , 

TeP 

T k :={TeP: \Z T \ = k} . 

The pairing CF for Z = m can be obtained by a variation 
of the procedure depicted above: 

lim (z 1 -2ei)(z 2 -2ei)C{£,<D,F) . 

Zi->2ej 

Z 2 -f2Ei 

But in the present case (sj = — 1/2 Vj) it is more con- 
venient employing the formula (|30| because S^Sj = 
1/2 ±5;. 

We comment that practical use of the formulas is lim- 
ited by the vastly increasing number of partitions, which 
depends on the number of pairs \£ \ = N and the order 
of the CF \H\. We want to emphasize that no complete 
knowledge of all the eigenstates is required. It doesn't 
show any dependence of the Hilbcrt space dimension ei- 
ther. We finally point out that the results apply to arbi- 
trary Sj (i.e. any degeneracy of the single particle levels). 



IV. DETERMINANT REPRESENTATION OF 
THE CORRELATION FUNCTIONS 

In this section wc will sketch how the formula obtained 
above can be simplified by recasting the CF's into sums 
of determinants of certain N x N matrices. As it was 
remarked above the scalar products between the BCS 
Bethe states (£\F) = C(£,0,J r ) can be expressed as de- 
terminants. Within the formalism we exploited in the 
previous section (see Ea.(|28"j)) in fact 



C{£,$,T) = {-1 



, A' 



E( \{n T (\£ T \)^S{£ T UF T ) \ , (34) 
where S(£ U F) can be written as 

Sklyanin proved]!?) that C(£,0,.F) can be indeed 
written as a polynomial that is linear in each of 
s(e) - s(f) 

S(e, f) = Ml e e £, f e J 7 (see Eq. ETJ. Con- 

e- / 

sistcntly with Richardson's old result |22j, it can be ex- 
pressed as a sum of TV! determinants 



C(£,0,F) = Y, M„ 



(36) 



where it is an element of the symmetric group Sn and 
Mir is defined as 

1 



(e a - ep)(fn(a) ~ fir(P)) 

A major simplification was achieved by Slavnov who was 
able to express the scalar product as a single determinant 
[29| . Therefore the Eq. ([36"]) can be recast into 



(£\7) = 



N N 

n u(u 

p=\ ai tp 



II (ep-e*) II (fp-U) 

P<a a<P 

det N H({f a },{ep}), (37) 

As discussed in SectHJ the CF's of the BCS model are 
identical to those of the Gaudin model, with the para- 
meters e , / satisfying Richardson's equations ([2"T]) in- 
stead of the Gaudin-equations. The entries of the N x N 
matrix H({f a }, {e^}) are 



H„ 



h 



E 



1 



1 



-2E 



if a - fa)(e b - /„) 



(38) 



6 



The norms of the states are obtained for e a f a in 
(|37| and give |V>(e Q )| 2 = deti^G^ where G is the Gaudin 
matrix given by 



Gab = { 2 



E 



; b = a 



(39) 



The various stages of the calculation of CF's proceed 
through certain recurrence formulas involving (|37|) as a 
basic ingredient (see [2l| for the details). Therefore the 
CF's result to be determinants as well. We comment that 
such a simplification was first achieved after a tour de 
force on integrable spin 1/2 theories (beyond the quasi- 
classical expansion) leading to the so called solution of 
the inverse problem [30l l3l| . The main accomplishment 
is that the lattice spin variables arc expressed in terms 
of the entries of the monodromy matrix A(u|e), B(u\e), 
C(u\e), D(u\e) in a closed (and simple) form. This al- 
lows to evaluate the CF's in the non-local Hilbert space 
spanned by the Bethe vectors (instead of expressing the 
Bethe states in Hi oc )- Zhou et al. [26j calculated the 
relevant quantities for the BCS model through the quasi- 
classical limit (fl2|). generalizing the solution of the in- 
verse problem to non-fundamental integrable spin theo- 
ries (where the auxiliary and the quantum spaces have 
different dimensions). The formulae read (2(| 



a— 1 a— 1 

i— 1 i 

a— 1 a— 1 

St = f[ t(e a )K-^ {A{et) - D{tl)) K^ t-\e a ) 



and Q given by 

N+l 



T ab (m) = J| (e„ - /&) E 



"3 1 \r- 

a^ta 



{ (fb ~ ej)(e - £j) 



-2E 



1 



(fb - e a )(e a - e a ) 



b<N + l, 



7aiv+i(m)= — — -, Qab(m) - 



(e a - e m ) 2 



(^a ^m) 2 



T is the N x N matrix obtained from T by deleting the 
last row and column and replacing TV + 1 by N in the 
matrix elements. Here, we assume that both {/<*} and 
{eb} are solutions to Richardson's Bethe equations (f2Tj) . 
The two-point correlation functions are 

{£\S-S+\F) = 
N 1 

E T - ^ 6l ' " ' ' eN \Sm\h> ' • ' J /«> • • • 7 /at) - 

Q = l /Q C ™ 



E 



1 



(f a ~ e n )(/8 - e„) 



(d, • • • ,e N \S m S n |/i, • • • ,/«,-•• , /a, • ' ' , In)- (42) 

Here, the hat denotes that the corresponding parameter 
is not present in the set. Since {e a } is a solution of the 
Bethe equations, (£|S^|/i, • • • , f a , • • • , /at) is the form 
factor given before, while 

n(£\S~S~ |J")jv-2 = 



N 

I! (e/3-e m )(e^-e„) - 
ff=i det N T{m 7 n, {epj, {,/«}) 

Uf A 11 (^ - e a )(/ T - / 4 ) 



(43) 



with 7\ := exp(— 2rj X^=i Sj/ffi) being essentially the 
total S z , and t(u) bein g th e transfer matrix. The form 
factors are obtained as KMJ 



Ar+i(f|5 m |J')iv - 
jN+i 



II/3=i ( e /3 - e m) detjy + iT(m, { e<9 }, {/„}) 

ria=l(/« - £ ™) ri/3>a(e/3 " O IT^a^ - fa) 



,(40) 



JV 



(d, ■ ■ • ,e N \S^\fi, ■■ • ,/jv) = J 



a=1 (/a — Cm) 



det w Uf ({e^}, {/„}) - Q(m, {e^}, {/ a }) 



n^> a ( e /3 - e «) ri/3< a (//3 - /«) 



. (41) 



where jy + i(£\ and |.F)jv indicate Bethe states with N + l 
and iV rapidities respectively. The matrix elements of T 



with 



N 



Tab(m,n) = Y[(e a f b ) ^ 



\ (fb - £j)(e - ej) 



i 



(/& - e Q )(e a - e Q ) 



, 6<JV-1, 



TaN-i(m,n)-- 
T a N(m,n)-- 



2c a 6 m e n 
[(ea - £m)(e a - e„)] 2 ' 
1 



(ea ^m) 2 



In (|43| m 7^ n is assumed, and it is zero if m = n. 

The expression for the charge CF can be obtained from 
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the form factors and CF above[32 



(£\S z m S* n \T) = 



(44) 



N 



+ 



5., 



iV 



2 2( e m/a) 2(e n / Q ) 



( e m — fa)( e n — f@) 



l/l: ' ' ' 7 /a, ' ' ' i /iv) + 

' • • Jo,"' j //Si ' ' ' 5 /at) 



The static CF of interest here are obtained with (£| = 



A. Reduction formulas 

The last significant progress for the evaluation of the 
CF for (£\ = (T\ was pursued by Faribault, Calabrese 
and Caux[32j . They managed to reduce the complexity of 
the above expressions to sums over only TV determinants. 

Both (S m S„) and (S^S*) involve the evaluation of 
the form factors (j4*0"j) and (|4"5j) for e — > /. In such a limit 
the CF are still a sum of two terms, each one involving 
sums of TV determinants of modified Gaudin matrices. 
In Ref. [32| the specific symmetry of the Richardson 
equations was exploited to reduce the CF to a single term 
expressed as a sum of TV determinants (see [HI for the 
detailed calculations) . In this way, the final result for the 
correlation function is 



N 



(£\S-S+\£) = J2 



_ D (m,n) ( 45 ) 



where 



D 



( 771,71 ) 

q,i 



Gi 



Ki, 



i+lq 



q=l 



■Gi+i 



Q + 2 (gg - fnKgg-1 - e m ) 



i<q-l, 
B i = q-l, 



Here, 



C 
Gi 



C a 

B n 



I = q, 
i > q. 



1 



(w a - £m) 2 ' 



(w a - e rn ) 2 (w a - £n) 2 

Similarly (S^S*) is given by [H 



(46) 

(47) 
(48) 



(£\s z m s*\£) 



N 



i£!£l - l^2(dctD ( ™' n) + &ctD^ m) ). (49) 

9=1 

The formulas obtained above for the correlation functions 
are completely general and are valid for any choice of the 



Hamiltonian parameters e Q and g (with some caveat in 
the limit of coinciding energies). The low level of com- 
plexity of this representation as sum of TV determinants 
of TV by TV matrices allows one to have access to the 
static correlation functions for systems with a reasonable 
number of pairs. 



V. CANONICAL PAIRING FLUCTUATIONS 

In the canonical ensemble the conventional BCS order 
parameter A = n(U\c^c^\U)n is vanishing exactly. Nev- 
ertheless the pairing instability can be characterized by 
studying the correlation function 

= (4Ai c H c n) - ( c lt c it)( c 4 c i4-) ( 50 ) 

indicating the tendency that electrons form Cooper pairs 
instead of uncorrelated electrons. The canonical BCS 
order parameter is[9l.l39l] 



n 



(51) 



which in the limit of large volume fl and large TV reduces 
to the BCS value (see Sect I VI [) . We observe that, in con- 
trast with a normal Fermi gas, a system with pairing 
instability will take an energetic advantage by increas- 
ing S7 for fixed SI /TV (because the phase space available 
for coherence is enlarged). Therefore, energy correlations 
are short ranged in a normal Fermi gas and long ranged 
in the presence of a pairing coherence. Accordingly, the 
footprint for an ongoing pairing instability is a finite size 
scaling ansatz 



(52) 



The above quantity was evaluated exploiting the ex- 
act formulas (|30|) and (|4Tj) . The rapidities involved in 
the equations at finite size are obtained by solving the 
Richardson equations numerically for the model parame- 
ters, which here are equally spaced single particle energy 
levels e l = i and half filling il = 2TV (see Ref. [Io[ for the 
details). The results for mVi are shown in Figs. [1] and 
[2j Whereas the former shows the g dependence at fixed 
TV, in the latter each plot consists of the various curves 
at fixed g (=0.1, 0.2, 0.4, 0.7) for varying TV. It is clear 
that the results tend to the BCS result. In Ref. [13] was 
noticed that this convergence is the slower, the smaller g 
is: for g = 0.1 the maximum at TV = 256 is only 90% 
close to the asymptotic result wheras at g = 0.7 the 
TV = 16 result is already at 99.8%. In Fig. [3] the or- 
der parameter vff /fl as a function of g for several values 
of TV is shown and compared with the BCS result: for 
TV = 128, ^ is almost indistinguishable from its limit- 
ing value for large enough g. The scaling of $ was 
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FIG. 1. (Color online) Pair distribution mVj, as function of i. 
Each plot is for fixed SI but for several different g going from 
0.1 to 1 in steps of 0.1. With kind permission by the authors 
of Hi. 
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FIG. 2. (Color online) Pair distribution UiVi as function of 
i/fl for different values of 0. The four plots are taken at 
variable g. With kind permission by the authors of [32T ]. 




FIG. 3. (Color online) Canonical order parameter ^ as a func- 
tion of g for different numbers of pairs S7 = N/2. The BCS 
result it = S7/(4g sinh(l/2(?)) is plotted for comparison. In- 
set: Small g behavior of ^ compared with large N expansion. 
With kind permission by authors of [32| . 
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FIG. 4. Scaling of Data for "J 7 , as a function of a = g/8 
and S7 = 2N — even, which show clearly a phase transition 
(N = 8(circles), N — 12(squares) , N = 16(diamonds) , N — 
20(triangles) ) . In the inset it is shown the data collapse. The 
parameters are gl = 0.315, n(gl) = 0.940, and \/v = 0.26. 
Taken from M- 



originally obtained in Q for small size at a given value 
of the pairing coupling g\ (see the caption of Figj4j) . In 
Ref . [28| a further scaling point gi was evidenced (see the 
caption of Fig{5j) . In order to extract the finite-size scal- 
ing, log[$n(g)/^Q>(g)]/ \og[il/Q'] was taken in consider- 
ation for different valus for Si and Si'. At a scaling point 
all these curves cross (r?(Sl, Si', g*) = r)(g*)) as shown in 
Fig. [5] The physical meaning of two apparent "scaling 
points" was unclear and deserved further analysis. This 



analyis was significantly extended in |32|, but without 
any scaling analysis. However, there is no second scal- 
ing point visible in their analysis for larger pair number 
(N > 32). In figure [6] we present the data collapse of the 
data of Ref. [32| . The value for r\ is in accordance with 
rj = 1 with an error about 0.01. The second coefficient 
v = 16.6 leads to the searched-for data collapse. 

It is instructive that the computed leading term of VP 
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FIG. 5. The figure shows two crossing points: the first is in 
agreement with Ref. [jj; the second is at g* — 0.417, g — 
1.028. The data collapse is seen in the inset for 1/v = 0.15. 
Taken from [H. 




FIG. 6. (Color online) The data collapse for rj = 1.00 ±0.001 
and v = 16.6. The data points have been read out of the 
figure 8 taken from Ref. [32 ]. 



for large f2 and g 

# _ 1 

N ~ 2 



+ 0(g-'\N- 1 ) 



(53) 



coincides with that of Ref. [36|. Instead, for small g and 
large TV 

| =g h(8W§)) +Q(1/]njy)[ forg<<1 _ (54) 



The N independent result for large g and the N^ 1 / 2 
dependence at small g gives a hint towards the non- 
pcrturbative nature of superconductivity. Further stud- 



ies on the finite size corrections of the BCS pairing am- 
plitude were performed in[41|. 

We close the discussion on the canonical pairing fluc- 
tuations mentioning the relation between ^ and the 
Onsager-Penrose-Yang parameter [42j for the long-range 
off-diagonal order^oD, taking into account the effect of 
non-diagonal correlations. 

Although $o£i can be obtained with the formulas for 
(S^~S~), a much easier route is to apply the Hellmann- 
Feynman theorem: 



*od = jz £ (S+Sr) = -± d E (g)/dg , 

i,j=l 



(55) 



where Eq is the ground state energy of the BCS model. 

In the thermodynamic limit ^qd 00 = ^"n °° ■ However, 
the two quantities are independent for finite sizes. Tian 
et al. [43| proved that \P and $ofl satisfy the following 
relations for any value of g and N 

- 1) < < 1 + . ( 56 ) 

For N — > oo these are trivial bounds, but not so for finite 
N. 



VI. THERMODYNAMIC LIMIT 

The Richardson equations (|2T|) admit an electrostatic 
analogy [l6|, 0, EEJ , where the eigenenergies £j and 
solutions e a both are interpreted as point charges of the 
strengths —1/2 and 1 respectively. The thermodynamic 
limit is performed making use of this analogy. 
Define 



p(Xj) 

a(z a ) 



1/2 



i 



^|e a+ i - e a \ 



(57) 
(58) 
(59) 



This choice leaves the Debeye-shcll invariant. Inserting 
this into the Richardson equation ([21]) , il cancels and we 
obtain 



1 

9 



ds 



\dz 



. 



(60) 



Following the works [l6|,ll8( the BCS-gap A BC s = 2 A is 
the imaginary opening of the arc solving the Richardson 
equations, a = A±iA (see Fig. [7]), and the gap equation 
becomes 



E 



de 



V(e 4 -A)2+A2 
d(e)p(e) 



V(e-A) 2 +A2 ' 



(61) 
(62) 
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FIG. 7. (Color online) Left: Solutions to the Richardson eqa- 
tion l21l Right: A picture on how the BCS-gap is determined. 



where di (d(e)) is the multiplicity of the level £j. This re- 
suls in the known expression A = sin h jy 2 g ^ or ^ ne ec l uaii y 
spaced model. A good summary of Gaudin's article to- 
gether with a modern numerical analysis of the solutions 
for the BCS-model is found in [Hj]. The generalization 
of the electrostatic analogy to more general settings (the 
tri gon ometric and hyperbolic BCS-model) can be seen 
in I45il. 



VII. FURTHER DIRECTIONS 

In this article we have reviewed the current under- 
standings of the mesoscopic fluctuations of the pairing 
instability based on Bethc ansatz techniques. The rele- 
vant quantity is a correlation function (CF), where the 
physical observables are evaluated as a static expectation 
value in the eigenstates of the Hamiltonian. The progress 
in the field of exact solutions are mature enough to allow 
an exhaustive analyis of superconductivity from meso- 
scopic regimes to thermodynamic limit at equilibrium. 
An important piece of information, however, comes from 
the study of the system out of equilibrium. The typical 
picture is provided by transport experiments were the 
dynamical CF, G(t) = (O(0)O(t)), are the interesting 
quantities to be calculated. The formula for G(t) involve 
an additional level of complexity. The basic ingredients 
are off diagonal correlations, namely static CF between 
different eigenstates. The first exact off-diagonal CF for 
the BCS-model obtained in [28| could only be calcu- 
lated for very small sizes. The better performance of 
the determinant expression in (2(| could be even fur- 
ther improved in |32| to make reasonably higher pair 



numbers accessible. Nevertheless, this is not the end of 
the story because in principle all the off-diagonal CFs 
are involved in G(t). Fortunately, the problem can be 
simplified for the BCS model because the eigenstates do 
not contain the coupling constant explicitly. In a very 
relevant paper, Faribault, Calabrese and Caux combined 
numerical and analytical analysis to realize that indeed 
only a relatively small amount of excitations contribute 
significantly to the dynamical CF (46|. As a consistency 
check they used exact sum rules relating the dynamical 
to static CF (the latter can be accessed easily) . They dis- 
covered that the weight of the multi-particle excitations is 
suppressed increasing N: in the thermodynamic limit the 
two-particle excitations are hence dominant in the calcu- 
lation of G(t) [47j. This is the ultimate reason why the 
Bogoliubov mean-field results (just neglecting the higher 
order correlations) coincide with the exact ones in the 
thermodynamic limit. 

An important problem that has been intensely studied 
in the recent literature is the response of a given sys- 
tem when pushed out of equilibrium by a sudden change 
in some control parameter: the quantum quench (see 
pH ] for a review). The richness and complexity of this 
problem is very much related to the developing nonlo- 
cal character of the correlations in the system by time 
evolution [49[. Remarkably, this kind of issues can be 
explored experimentally at the quantum level by realiz- 
ing highly controllable quantum many-particle systems 
with cold atoms (Hoj . We observe, however, that pairing 
fluctuations in cold atoms are expected to be more evi- 
dent in transport experiments rather than in the popular 
expansion protocols (where the increase of single parti- 
cle kinetic energy might mask the crossover). Arrays of 
coupled microcavities are potentially interesting alterna- 
tive experimental platforms [HI, (HJ • Those problems are 
studied through the dynamical CF as well. The time 
evolution starts, because the eigen-basis where the wave 
function of the system lives, changes after the quantum 
quench. The computational complexity of the problem, 
generically, increases factorially with the size of the sys- 
tem. The problem of the quench dynamics in the in- 
tegrable BCS model, (Eq. ([!]), was thoroughly studied 
in (HI]; see also [HI). By employing the approach de- 
veloped in [46[ the authors proved, first, that the all the 
quench matrix is accessible by the Slavnov formula; then 
they proved that the quench dynamics occurs only along 
a relatively small subspace in the Hilbert space. Athough 
their results provide a hint that deviation from the mean- 
field regime emerges in the quench dynamics at finite size, 
further analysis seems to be required to unambiguously 
disclose the effect of mesoscopic pairing fluctuations. 
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